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The first generation of stars had very different properties than later stellar genera- 
tions, as they formed from a "pristine" gas that was completely free of heavy elements. 
Normal star formation took place only after the first stars polluted the surrounding tur- 
bulent interstellar gas, increasing its local heavy element mass concentration, Z, beyond 
Q-f a "critical" threshold value, Z c (10 -8 Z c <> 10~ 5 ). Motivated by this astrophysical 

problem, we investigate the fundamental physics of the pollution of pristine fluid ele- 
ments in statistically homogeneous and isotropic compressible turbulence. Turbulence 
stretches the pollutants, produces concentration structures at small scales, and brings 
the pollutants and the unpolluted flow in closer contact. The pristine material is pol- 
luted when exposed to the pollutant sources or the fluid elements polluted by previous 
mixing events. Our theoretical approach employs the probability distribution function 
(PDF) method for turbulent mixing, as the fraction of pristine mass corresponds to the 
low tail of the density-weighted concentration PDF. We adopt a number of PDF closure 
models and derive evolution equations for the pristine fraction from the models. To test 
and constrain the prediction of theoretical models, we conduct numerical simulations for 
decaying passive scalars in isothermal turbulent flows with Mach numbers of 0.9 and 6.2, 
and compute the mass fraction, P(Z C , t), of the flow with Z ^ Z c . In the Mach 0.9 flow, 
the evolution of P(Z c ,t) is well described by a continuous convolution model and goes 
as P(Z c ,t) = P(Z c ,t)ln[P(Z c ,t)]/T con , if the mass fraction of the polluted flow is larger 
than ss 0.1. If the initial pollutant fraction is smaller than w 0.1, an early phase exists 
during which the pristine fraction follows an equation derived from a nonlinear integral 
model: P(Z c ,t) — P(Z c ,t)[P(Z c ,t) — l]/r int . The timescales r con and r int are measured 
from our simulations. When normalized to the flow dynamical time, the decay of P(Z C , t) 
in the Mach 6.2 flow is slower than at Mach 0.9 because the timescale for scalar variance 
decay is slightly larger and the low tail of the concentration PDF broadens with increas- 
ing Mach number. We show that P(Z c ,t) in the Mach 6.2 flow can be well fit using a 
formula from a generalized version of the self-convolution model. 



1. Introduction 

Big bang nucleosynthesis produced helium efficiently, but it was halted by the expan- 
sion of the universe before it was able to make stable elements heavier than lithium 
(Walker et al. 1991). On the other hand, even the most pristine stars observed today 
(Cayrel et al. 2004; Frebel et al. 2008; Caffau et al. 2011) have substantial mass fractions 
of heavier elements, indicating that they have been polluted with the nucleosynthesis 
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products of an as-yet undetected first generation of stars. This early stellar generation 
had an enormous impact on the evolution of later forming stars, and such stars are likely 
to have been much more massive (Abel, Bryan, & Norman 2000; Bromm, Coppi & Lar- 
son 2002) and much hotter (Schaerer 2002) than present-day stars, due to the important 
role that heavy elements play in star formation and evolution. 

When and where this remarkable early stellar generation formed is a question of fun- 
damental astrophysical importance. On cosmological scales, the key issue is the time it 
takes for heavy elements to propagate from one galaxy to another. As shown in Scan- 
napieco et al. (2003), the distances between these regions of early star-formation are so 
vast that the universe was divided into two regions: one in which galaxies formed out of 
material that was already polluted with heavy elements and one in which galaxies were 
formed from initially pristine material. 

This second set of initially-pristine galaxies is especially interesting, as the first stars 
formed in these galaxies may be observable (Scannapieco et al. 2005; Jimenez & Haiman 
2006; Nagao et al. 2008). As star formation continued in these objects, the interstellar gas 
became enriched with heavy elements released by the explosions of the first stars. This 
self-enrichment process increased the abundance or mass fraction, Z, of heavy elements, 
and finally led to a transition to normal star formation in regions where Z exceeds a 
critical value, Z c . This critical value is expected to lie in the range 10~ 8 <, Z c <, 10~ 5 
(or 10~ 6 to 10~ 3 times the heavy-element abundance in the Sun), depending on whether 
the cooling of the interstellar gas is dominated by dust grains (Omukai et al. 2005) or by 
the fine structures of carbon and oxygen (Bromm & Loeb 2003). 

In a given galaxy, the key quantity to characterize the transition to normal star for- 
mation is the fraction, P(Z c ,t), of the interstellar gas with Z below Z c as a function of 
time. The temporal behavior of this fraction depends not only on the rate at which new 
sources of heavy elements are released to the interstellar gas, but, more importantly, on 
the transport and mixing process of these elements in the galaxy (Pan & Scalo 2007). For 
example, a high mixing efficiency would result in a rapid decrease in P{Z Cl t), and hence 
in a sharp transition as the average concentration of heavy elements exceeds the thresh- 
old Z c . On the other hand, a low mixing efficiency would lead to a gradual transition. 
The interstellar gas in these galaxies is expected to be turbulent and highly compress- 
ible, and the turbulent motions are likely to be supersonic (Grief et al. 2008; Wise et al. 
2008). Therefore, understanding mixing in supersonic turbulence is crucial to answering 
the question of how the pristine gas in early galaxies was polluted. 

In the present paper, we do not intend to directly model the complicated mixing pro- 
cess in a realistic galactic environment. Instead, we investigate the fundamental physics 
of turbulent mixing in compressible flows using idealized analytical and numerical tools. 
The primary goal is to understand the pollution of pristine material in statistically ho- 
mogeneous and isotropic turbulence. This underlying physics is prerequisite for modeling 
the mixing of primordial gas in realistic interstellar turbulence. In a future work, we will 
apply the results of the current study to build a subgrid model for large-scale simulations 
for the formation and evolution of early galaxies. These simulations account for the com- 
plexities in the interstellar medium, but cannot resolve the scales at which true mixing 
occurs. The subgrid model will provide a crucial step toward predicting the transition 
from primordial to normal star formation in the first generation of galaxies. 

A systematic numerical study of passive scalar physics in supersonic turbulence has 
been recently conducted by Pan and Scannapieco (2010, 2011), who simulated scalar 
evolution in six compressible turbulent flows with Mach number ranging from 1 to 6. In 
these papers, a detailed analysis of various statistical measures for the scalar field was 
performed, including the scalar dissipation, the scalar probability distribution, the power 
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spectrum, the structure functions and intermittency. It was found that the classic cascade 
picture for passive scalars in incompressible turbulence is generally valid also for mixing 
in supersonic turbulent flows. The effect of compressible modes in supersonic turbulence 
and their modifications to the classic picture for passive scalar turbulence were examined 
by analyzing the Mach number dependence of the scalar statistics. The conclusions of 
these studies provide general theoretical guidelines for understanding the mixing process 
in interstellar turbulence. 

To explore how the pollutant-free mass is contaminated in turbulent flows, we make use 
of the probability distribution method for turbulent mixing. The fraction of unpolluted or 
slightly polluted flow mass corresponds to the far left tail of the probability distribution 
function (PDF) of the concentration field, as Z c is typically much smaller than the average 
value. This fraction can be evaluated by integrating the concentration PDF from zero to 
the threshold, Z c . We will generally refer to the fraction P(Z C , t) as the pristine fraction. 
Note that our approach here is general, and is not limited to mixing in early galaxies. 

The PDF equation for passive scalars cannot be solved exactly because of the closure 
problem, and various closure approximations have been developed to model the PDF 
evolution. In this work, we consider several existing closure models and derive equations 
for the fraction P(Z C , t) for each of them. The far left PDF tail corresponds to high-order 
moments of the PDF, and thus it is quite uncertain whether the closure models can 
capture the high-order statistics with sufficient accuracy. In order to test the reliability 
of the adopted models and constrain their parameters, we perform numerical simulations 
for turbulent mixing in a transonic flow and a highly supersonic flow. 

The structure of this paper is as follows. In §2, we present the general PDF formulation 
for mixing in compressible turbulence. §3 gives a brief description for several existing 
closure models for the diffusivity term in the PDF equation. The predictions of these 
models for the mass fraction of unpolluted or slightly polluted flow are derived in §4. 
We describe our numerical simulations in §5, which are used to test and constrain the 
theoretical models in §6. Our main conclusions are summarized in §7. 

2. PDF Formulation for Mixing in Compressible Flows 

The PDF formulation was first developed for the probability distribution of the tur- 
bulent velocity field by Monin (1967) and Lundgren (1967), and for the PDF of the flow 
vorticity by Novikov (1967). The derivation of Monin (1967) was based on the equation 
for characteristic functions of the velocity field, while Lundgren (1967) started directly 
from the conservation laws of the flow. The two methods were later extended to derive 
PDF equations for scalar fields convected by a turbulent flow, such as the flow temper- 
ature or enthalpy, and the concentration fields of passive or reactive species in the flow 
(Ievlev 1973; Dopazo and O'Brien 1974; Pope 1976; O'Brien 1980; Pope 1985; Kolle- 
mann 1990; Dopazo ct al. 1997). Recent discussions of PDF equations for passive or 
active scalar turbulence can be found in the monograph by Fox (2003) and the thorough 
reviews by Veynante & Vervisch (2002) and Haworth (2010). 

In Appendix A, we derive the PDF equation for a passive scalar in compressible tur- 
bulence using the method of Lundgren (1967). The derivation is based on the equation 
of the concentration field, C(x, £), for passive tracers advected in a turbulent flow with 
density p(x,t) and velocity v(x, t): 



where the concentration field is defined as the ratio of the local tracer density to the 



dC_ 

dt 



+ v • VC* = -V • (pnVC) + S*(x, t), 



P 



(2.1) 
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flow density. In the diffusion term, k denotes the kinematic molecular diffusivity, and the 
dynamic diffusivity, pn, is basically independent of p (i.e., k oc p^ 1 ). The term S^x, t) 
represents the sources of new pollutants. 

Our derivation in Appendix A adopts a density-weighting scheme, which is appropriate 
for passive scalar mixing in compressible flows (Pan and Scannapieco 2010). We define a 
density-weighted concentration PDF, p(Z; x, t) = (p5[Z—C(x, t)]), where (•••) denotes the 
ensemble average, the density- weighting factor p = p(x.,t)/p is the ratio of the local flow 
density to the average density p, and Z is the sampling variable. Using the advection- 
diffusion equation (|2.ip . and the continuity equation for the evolution of the density- 
weighting factor, we obtain, 

0p(Z;x,t) / (pv\C = Z) \ = d ( (V ■{ P kVC)\C = Z) \ d ( ( P S\C = Z) 
dt + \ P {p\C = Z)J dz\ P (p\C=Z) ) dz\ P {p\C = Z) 

(2.2) 

where {■ ■ - \C = Z) denotes the ensemble average under the condition that the concentra- 
tion field C(x, t) is equal to Z (see Appendix A). The equation is essentially a Liouville 
equation for the conservation of the concentration probability. To our knowledge, this 
equation for the scalar PDF with density-weighting has not been derived before. 

The density weighting scheme is preferred in our study for two reasons. First, rather 
than the volume fraction, we are interested in the mass fraction of pristine gas in early 
galaxies, which corresponds to the left tail of the density-weighted PDF. Second, the ad- 
vection term in the equation for the density- weighted PDF takes the form of a divergence, 
and thus conserves the global PDF (i.e., the integral of the local PDF, p(Z;x, t), over 
the flow domain). This provides a formal and rigorous proof for the physical intuition 
that the turbulent velocity field itself does not homogenize the distribution of pollutants. 
The advecting velocity transports, redistributes and deforms the concentration field, but 
does not change the mass fraction of fluid elements with a given concentration level. 
Furthermore, the advection term vanishes if the flow and the concentration fluctuations 
are statistically homogeneous. 

In contrast, if one derives an equation for the volume-weighted PDF for a passive 
scalar in compressible turbulence, the advection term would not be a divergence term. 
The term reflects the effect of flow compressions and expansions, which can change the 
volume fraction of fluid elements at a given concentration (Pan & Scannapieco 2010). 
This effect on the PDF is clearly different from scalar homogenization, and can be avoided 
by adopting a density-weighting factor. We thus argue that it is more appropriate to use 
the density-weighted PDF equation for the study of mixing in compressible turbulence. 

Molecular diffusion is the only process that homogenizes, and the molecular diffusiv- 
ity term in the PDF equation continuously reduces the PDF width. This term can be 
rewritten as, 



d_ ( (V • {pnVC)\C = Z) 
dz\ P ( P \C = Z) 



d_ t { P nVC\C = Z) 
dz\ p ( P \C = Z) 



d 2 ( { P k(vc) 2 \c = zy 
p- 



dZ 2 V (p\C = Z) 
(2.3) 

where both terms on the r.h.s depend on the ensemble average of the concentration 
gradients conditioned on C(x, t) — Z. As it is a divergence term, the first term in eq. 
(|2.3p conserves the global concentration PDF. The scalar homogenization is achieved 
through the second term, which is essentially a diffusion term with a negative coefficient in 
concentration space. This term keeps narrowing the concentration PDF, and the physics 
of turbulent mixing can be viewed as an anti-diffusion process in concentration space. 
Taking the second order moment of the last term in eq. (|2.3|) gives the scalar dissipation 
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rate, — 2(/5k(VC) 2 ). Using this rate, we define a mixing timescale, 

(~p5C 2 ) 
Tm ~ 2( / 5k(VC) 2 )' 



(2.4) 



where SC = C — (pC) is the fluctuating part of the concentration field. The timescale, 
T m , corresponds to the scalar variance decay by mixing, and thus characterizes the rate 
at which the diffusivity term reduces the PDF width. Although the diffusivity terms 
in eqs. (|2.3p and (|2.4I) do not have an explicit dependence on the flow velocity, the 
mixing timescale is determined primarily by the turbulent velocity field. This is because 
the turbulent velocity produces progressively smaller structures and thus strongly ampli- 
fies the scalar gradients, (VC) 2 , in eqs. (|2 .3[) and (|2.4p . By feeding molecular diffusivity 
with large-gradient structures, turbulent motions greatly accelerate the scalar dissipa- 
tion/homogenization. 

In the classic phenomenology for mixing in incompressible turbulence, the generation 
of small-scale concentration structures is through a cascade process similar to that of 
kinetic energy (Obukohov 1949; Corrsin 1951). The cascade is caused by continuous 
turbulent stretching, and it starts from the scale where the pollutant sources are injected 
into the flow, and proceeds to the diffusion scale where the molecular diffusion efficiently 
homogenizes the scalar fluctuations. The diffusion scale is essentially the scale where 
the action of molecular diffusivity becomes faster than turbulent stretching. From this 
picture, the mixing timescale, r m , is determined by the cascade time, which is essentially 
the eddy turnover time at the injection scale of the scalar sources because the cascade 
becomes faster and faster with decreasing length scale. 

Pan & Scannapieco (2010) showed that the cascade picture also applies for mixing 
in supersonic turbulence. They found that the mixing timescale, r m , was close to the 
eddy turnover time at the pollutant injection scale in all their simulated flows with Mach 
numbers in the range from 1 to 6. The existence of compressible modes in supersonic 
flows causes only a slight Mach number dependence of the mixing timescale, and the 
primary "mixer" is the stretching by solenoidal modes even at very high Mach numbers. 

Translating the physical discussion above to the mixing process of the unpolluted fluid 
elements in a turbulent flow gives the following picture. The turbulent velocity stretches 
the pollutants into smaller and smaller structures, and brings them to closer contact 
with the unpolluted flow. When the separation between the pollutant structures and the 
unpolluted fluid elements becomes close to or smaller than the diffusion scale, molecular 
diffusivity efficiently mixes them, reducing the unpolluted mass fraction. This suggests 
that the timescale for turbulent mixing to contaminate the unpolluted mass is also on 
the order of the scalar cascade timescale. 

The diffusivity term in the PDF equation has to be approximately modeled because 
of the closure problem (e.g., Dopazo and O'Brein 1974). Extensive efforts have been 
made to develop closure models for this term, and we will use several existing models 
in the current study, as described in §3. The advection term also has a closure problem, 
but modeling this term is not necessary if the flow and the scalar field are statistically 
homogeneous. The last term in the PDF equation ()2.2[) corresponds to the effect of the 
pollutant sources. In reacting turbulent flows, the source term due to chemical reactions 
has a closed form in the PDF formulation (e.g., Pope 1976), and this has led to the 
wide use of the PDF method in studies of chemical reactions in turbulent flows. In the 
present study, the pollutant source is merely an initial scalar condition, and we do not 
discuss modeling the source term further (see Pan and Scalo 2007 for an example with a 
persistent source). 
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3. PDF Modeling 

3.1. General Approach 

We employ both theoretical and numerical tools in the present work. The PDF formu- 
lation in §2 provides a general theoretical framework, and, to complete the theoretical 
approach, we will consider several existing closure models for the diffusivity term in eq. 
(|2.2j) . We will compare the predictions of these models for the scalar PDF evolution with 
that measured from numerical simulations. From our simulation data, we compute the 
concentration PDF as p(Z;t) = y J v p5[Z — C(x, t)]dx, where V is the total volume of 
the simulation box. This PDF measures the concentration fluctuations over the entire 
flow domain, and thus should be viewed as a global PDF. As pointed out in §2, the advec- 
tion term conserves the global density-weighted PDF, and thus need not be considered in 
our tests of the theoretical models against simulations. The global PDF is expected to be 
equal to the local PDF, p(Z; x, t), defined in the ensemble context under the assumption 
of statistical homogeneity. 

In our simulations, we only evolve decaying scalars with the source term 5(x, t) set to 
be zero. Neglecting the advection and source terms, the PDF equation becomes 

dp(Z;t) _ 3 f p (V-( P KVC)\C = Z) \ 



dt dZ V (p\C = Z) 

The only term that contributes to the PDF evolution in our simulations is the diffusivity 
term, and modeling this term is the main task of the PDF approach to turbulent mixing. 
In incompressible turbulence, the flow density is constant, and eq. (|3.1[) reduces to, 

^§^ = -^(p(^C)\C^Z)), (3.2) 

which has been extensively studied and modeled. 

The second order moment of eq. (|3.1I) corresponds to the scalar variance equation, 

al(5Z 2 ) _ (5Z 2 ) 



dt 



(3.3) 



where (SZ 2 ) = (pSC 2 ) denotes the density- weighted variance, and we have used the 
definition, eq. (|2.4[) . of the mixing timescale. In terms of the PDF, the variance is given 
by (SZ 2 ) = J{Z- Z) 2 p(Z, t)dZ with Z = J Zp(Z; t)dZ being the mean concentration. In 
general, r m may be a function of time. But if it is constant, the scalar variance decreases 
exponentially, which is the case at the late evolution stage of a decaying scalar (see §6.3). 

In analogy to the enrichment of pristine gas by the first generation of stars, the initial 
condition of the decaying scalars in this study will be set to be bimodal: consisting of 
pure pollutants (Z — 1) and completely unpolluted flow (Z = 0). This corresponds to a 
double delta function form for the initial concentration PDF, 

p(Z;0) = P 5(Z) + P 1 5(Z-l), (3.4) 

where P\ and Pq are the initial mass fractions of the pollutants and the unpolluted flow, 
respectively, and we have Po + Pi = 1 from the normalization of the PDF. 

The rest of this section is devoted to modeling the diffusivity term in the PDF equation. 
A variety of closure models have been proposed for this term, and the interested reader 
is referred to Dopazo et al. (1997) and Haworth (2010) for reviews. Here we will consider 
three of the models proposed for the diffusivity closure: the mapping closure model by 
Chen at al. (1989), the nonlinear integral models by Curl (1963), Dopazo (1979) and 
Janicka et al. (1979), and the self-convolution models by Villermaux & Duplat (2003), 
Venaille and Sommeria (2007) and Duplat & Villermaux (2008). 
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We point out that, in compressible turbulence, the diffusivity term has an explicit 
dependence on the density field, or more precisely, on the joint statistics of the density 
and the concentration fields. Therefore, an ideal PDF model for mixing in supersonic 
turbulence needs to account for the effect of density fluctuations on the diffusivity term, 
and to predict the dependence of the concentration PDF on the flow compressibility. 
However, to our knowledge, this has not been considered in existing models, which were 
usually tested against simulation results for mixing in incompressible turbulence. We will 
compare the predictions of the closure models mentioned above with our simulation data, 
and examine whether, by adjusting their parameters, these models can be successfully 
applied to the contamination process of pollutant-free mass in compressible turbulent 
flows at different Mach numbers. Future studies are motivated to develop closure models 
that explicitly address the effects of shocks and the Mach number dependence of the 
passive scalar PDF in supersonic turbulence. 

3.2. The Mapping Closure Model 

We first discuss the mapping closure model developed by Chen et al. (1989) for mixing 
in incompressible turbulence. We give a brief introduction to the model, and a detailed 
derivation can be found in, e.g., Pope (1991). The model is based on a surrogate field, 
(f>(x, t), obtained from the mapping of a Gaussian reference field #(x, i), 

^(x,i) =X[0(x,i),i], (3.5) 

where X is an ordinary, non-stochastic function. The Gaussian field, 0(x, t), is assumed to 
be statistically homogeneous and have zero mean and unit variance, i.e., the probability 
of finding 6*(x, t) equal to a given value n is given by g(rj) = — i== exp(— rj 2 /2) (see Girimaji 
(1992) for a generalized version of the mapping closure where the reference field PDF is 
time-evolving and not limited to Gaussian). The main idea of the model is to pursue a 
mapping function, X(r],t), with which the PDF of the surrogate field obeys exactly the 
same equation [i.e., eq. (|3.2[) ] as the actual field, C(x,t). This is indeed achieved if the 
mapping function evolves as 

dX(rj,t) k fd 2 X( V ,t) dX(r,,t) \ 

v — 5 — < ( 3 - 6 ) 



dt X 2 e {t) \ dif ' dn 

where X$(t) = ((V0) 2 ) -1 ' 2 , and Xg(t)/n is a timescale that controls the rate at which 
the mapping function and hence the PDF evolve. The timescale is unspecified in the 
original model, and can be calibrated by a comparison of the variance decay in the model 
with simulation results (see He and Zhang (2004) for a theoretical evaluation of this 
timescale using a two-point closure strategy). The evolution equation for the mapping 
function was derived in the incompressible limit, and thus the model is intended for 
mixing in incompressible turbulence only. By a comparison with our simulation data, we 
will examine whether the mapping closure model may also give acceptable predictions 
for mixing in compressible turbulence. 

With the desired mapping function, one can approximate the PDF of the actual field 
with that of the surrogate field. Using eq. (|3.6p . the PDF of the surrogate field can be 
converted from the Gaussian PDF of the reference field. The conversion gives, 

P(z;t) = 9iv )(^-y\ (3.7) 

where rj is the solution of X(r),t) = Z. 

The linear equation for the mapping function, eq. (|3.6[) . can be solved analytically, 
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provided the initial condition X(rj,0) (Pope 1991). For a double-delta initial PDF, the 
initial mapping is a Heaviside step function X(r], 0) = H(i] — 770), where rjo satisfies 
I-°oo 9{ 7 l)dri = Pq. With this initial condition, X(rj, t) is solved by 



x(n,t) = G 



E(t) 



(3.8) 



where T,(t) 2 = exp[JL K/\g(t')dt'] — 1, and G(rj) = f2 oo g(v')dri' is the cumulative func- 
tion of the Gaussian function. Combining eqs. (|3.7I) and (I3.8|) gives the predicted PDF 
evolution by the mapping closure model. 



3.3. The Nonlinear Integral Models 

In this subsection, we consider a class of closure models that use an integral form to 
approximate the diffusivity term in the PDF equation. This type of models originates 
from the equation introduced by Curl (1963), 



dt 



lit) 



dZ lP (Z i; t) / dZ 2 p{Z 2 ;t)5 Z 



p(Z;t)}, (3.9) 



where j(t) is the turbulent stretching rate. A physical interpretation of this equation is as 
follows (e.g., Dopazo 1979). The turbulent velocity field stretches the concentration field 
and produces structures at small scales. As shown by various studies, these structures 
are primarily in the form of 2D sheets (e.g., Pan and Scannapieco 2011 and references 
therein). The scalar sheets are brought closer to each other over time by the turbulent 
velocity. When the typical width and separation of the sheets decrease to the diffusion 
scale, molecular diffusivity can operate efficiently and homogenize. The timescale for this 
process is 7(f) -1 , which is expected to be on the order of the scalar cascade timescale or 
the mixing timescale r m . Two scalar sheets of different concentrations brought to close 
contact are assumed to mix perfectly, resulting in a concentration value equal to their 
average prior to the mixing event [see the delta function in eq. (|3.9[) ]. The last term in 
eq. (|3.9I) corresponds to the "destruction" of the previous PDF by the mixing event. 

One problem of Curl's model for turbulent mixing is that, if the initial concentration 
PDF consists of two delta functions [see eq. (|3.4j) ] , the predicted PDF shows unphysical 
spikes in between the initial delta functions. To avoid this problem, Dopazo (1979) and 
Janicka et al. (1979) independently generalized Curl's model replacing the delta function 
in eq. (|3.9[) by a smooth function J(Z; Zx, Z 2 ), 



d P (Z;t) 
dt 



lit) 



1 1 
p(Z i; t) J P {Z 2 ;t)J(Z;Z 1 ,Z 2 )dZ 1 dZ 2 




p(Z;t) 



(3.10) 



where JiZ; Z\, Z 2 ) represents the effect of mixing between two nearby scalar sheets with 
concentration values of Z\ and Z 2 . The function J(Z; Zx, Z 2 ) is zero for Z outside the 
range {Zx, Z 2 ) (or (Z 2 , Zx) if Zx > Z 2 ), and its normalization is J^ 2 J(Z; Zx, Z 2 )dZ = 1. 
A simple assumption for J{Z; Zx, Z 2 ) is that it is uniform between Zx and Z 2 , leading to 



d P jZ;t) 
dt 



= 7(*) 



z x 

p(Zx;t) 

z 



Zx 



p(Z 2 ;t)dZxdZ 2 



-p{Z-t) 



(3.11) 



where JiZ; Zx,Z<2) was set to 1/\Z 2 — Z\\ (Dopazo 1979 and Janicka et al. 1979). 
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The parameter j(t) as a function of time can be fixed by comparing the variance equa- 
tion of these models with the simulation data. The derivation for the variance equation 
can be found in Janicka et al. (1979) or Valino and Dopazo (1990). For Curl's model and 
the model with uniform J(Z; Zi, Z 2 ), the variance decays as oc exp[— | f Q -y(t')dt'] and 
exp[— -i J Q j(t')dt'], respectively. If the variance decreases exponentially with a constant 
timescale r m , y(t) is constant and equal to 2/r m and 3/r m , respectively, for the two 
models. 

Pope (1982) pointed out a weakness of this class of models: the normalized high-order 
moments, (8Z m ) / (SZ 2 )" 1 / 2 , do not converge with time for m ^ 4 (see also Valino and 
Dopazo 1990). This suggests that the predicted PDF by these models has excessively fat 
tails at late times (Kollemann 1990). 



3.4. The Self- convolution Models 

The last type of models we consider are those based on the self-convolution of the scalar 
PDF, which can be viewed as extensions of the model by Curl (1963) in Laplace space. A 
review for the development of these models can found in Duplat and Villermaux (2008). 

The Laplace transform t) of the scalar PDF is defined as p(£; t) = f Q p(Z;t)exp(—ZQdZ. 
Using the convolution theorem, the Laplace transform of eq. (|3.9I) gives, 

^M= 7 [j5(C/2;i) 2 -p(C;0]- (3-12) 

A similar equation in Fourier space was used by Pumir et al. (1991). This equation shows 
that turbulent mixing is essentially treated as a self-convolution process in Curl's model. 
We rewrite eq. fl3.12[) in a difference form p((; t + St) = ep((/2; t) 2 + (1 — e)p((; t) where 
e = 7<5i with St an infinitesimal time step. The difference equation can be interpreted as: 
during a time step St, mixing occurs only in an infinitesimal fraction, e, of the flow, and in 
this part of the flow the scalar PDF undergoes a complete convolution. The convolution 
process in eq. (|3. 12|) appears to be "discrete" . 

Following Venaille and Sommeria (2007), we derive a continuous version of Curl's 
model. We assume that, in each time step St, the PDF convolution occurs everywhere in 
the flow, but the number of convolutions is taken to be infinitesimal and equal to e (Duplat 
and Villermaux 2008). This assumption can be written as p((; t+St) = p(C/ (l + e ); t)( 1+€ \ 
Using the Taylor expansion p((/(l+e); i)( 1+e ) ~ p((; t)+e\p((; t) ln(p(C; t))-C,dp{£; t)/8Q 
and taking the limit St — > 0, we obtain, 



dp((;t) 
dt 



pln(p) -c|| 



(3.13) 



which represents the model of Venaille and Sommeria (2007). We will refer to this 
model as the continuous convolution model. In this model, the variance decays as oc 
exp(— J Q "f(t')dt r ). Venaille and Sommeria (2007) showed that the predicted PDF by eq. 
(|3.13p evolves toward Gaussian in the long time limit (in contrast to the integral models 
in §3.3). A comparison of this model with experimental data is given in Venaille and 
Sommeria (2008). We note that, if the initial PDF is two delta functions, the continuous 
self-convolution model is not applicable for the PDF evolution right from the beginning 
(Venaille and Sommeria 2007). We thus cannot compare the model prediction for p(Z; t) 
with our simulation results at the early evolution stage. The model will only be used to 
study the evolution of the unpolluted mass fraction. 

As pointed out by Duplat and Villermaux (2008), a more general extension of Curl's 
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model is, 

5p(C;<) 



at =in 



(3.14) 



Curl's original model (eq. (|3.12[1 ) and the model of Venaille and Sommeria (2007) (eq. 
()3 . 1 3[) ^) are special cases of eq. (|3.14[) with n = 1 and n — > oo, respectively. The parameter 
n can be a function of time in general. The assumption behind eq. (|3.14|) is that a fraction, 
ne, of the flow experiences mixing/convolution events during a time step St, and the 
number of convolutions in this fraction of the flow is 1/n. For eq. (|3.14[) . the variance 
decay goes like oc exp(— J y(t')n/(n+l)dt'). We will refer to eq. (|3.14[) as the generalized 
convolution model. 

We finally consider the model by Villermaux and Duplat (2003), which was motivated 
by a turbulent mixing picture with three related processes: the generation of pollutant 
sheets by turbulent stretching, the diffusion of the pollutant sheets by molecular diffu- 
sivity and the merging of the diffused sheets. The merging of the sheets corresponds to 
a self-convolution process. The model is represented by (Duplat and Villermaux 2008), 



dp(C,t) dp 

dn(t) 
dt 



p(C;*) (1+1/n) -i5(C;*) , 

= 7 n, (3.15) 



where n increases with time and the first equation can be viewed as the expansion 
of eq. (|3.14[) at large n. Note that eqs. (|3.15p and (|3. 131) approach the same limit as 
t — > oo. Villermaux and Duplat (2003) showed that eq. (|3.15l) has an asymptotic solution 
p(C;i) = (1 + (Z)n)~ n a t l ar S e t, which corresponds to a Gamma distribution for the 
scalar PDF (Duplat and Villermaux 2008), 

n n / nZ\ 

^^w^ Iexp U)' (3 ' 16) 

where T(n) is the Gamma function. The Gamma distribution is valid only at late times 
with n J> 1, and cannot be applied to study the pristine mass fraction at the early 
evolution stage when the fraction is significant. Therefore, we do not use the model for 
the pristine mass fraction, but will check whether the scalar PDF in our simulations 
approaches a Gamma distribution at late times. 

We point out a fundamental difference between the mapping closure model discussed 
in §3.2 and the models presented here and in §3.3. The mapping closure is established by 
a direct approximations of the exact, but unclosed form of the diffusivity term. On the 
other hand, the nonlinear integral models and the convolution models do not start from 
the diffusivity term in the PDF equation, instead they are largely based on a physical 
picture for the mixing process. 



4. Mass Fraction of Unpolluted or Slightly Polluted Flow 

As mentioned in the Introduction, we are interested in the mass fraction, P(Z c ,t), of 
the flow with concentration smaller than a tiny threshold, Z c , which can be calculated 
from the concentration PDF as 

P(Z c ,t) = J p{Z'-t)dZ'. (4.1) 
o 
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The fraction corresponds to the far left tail of the PDF since the threshold Z c of inter- 
est is typically much smaller than the average concentration. Taking the limit Z c — > 
in eq. (I4.1[) . we obtain the fraction P(t) of exactly pollutant-free mass, i.e., P{t) = 
limz c ->.o J Q c p(Z';t)dZ' . This fraction is zero unless p(Z;t) has a delta function compo- 
nent, S(Z), at Z = 0. In this section, we derive equations for P(Z c ,t) and P(t) from the 
closure models discussed in §3. 

An interesting observation of the action of molecular diffusivity is that it tends to de- 
crease the exactly pollutant- free fraction, P(t), to zero instantaneously. For illustration, 
we consider a simple situation with a point source diffusing in a static uniform medium. 
The concentration field obeys the diffusion equation, whose solution is given by a Gaus- 
sian function. From this solution, it is clear that, no matter how small the molecular 
diffusivity, k, is, the concentration field at a finite time (t > 0) becomes nonzero at any 
finite distance (r < oo) from the initial source, suggesting that all the pollutant-free mass 
is removed from the system instantaneously. 

This acausal behavior of molecular diffusivity originates from the Laplacian operator 
in the diffusion equation, which implicitly assumes that the random walk of some tracer 
molecules can bring them to an infinite distance during any small (but macroscopic) 
time interval. This is clearly unrealistic. The thermal motions of tracer molecules must 
have a finite maximum speed, max(w t h), and thus none of them can reach an infinite 
distance instantaneously. If the size of the system in question is L, there could be exactly 
pollutant-free mass surviving for a finite time ~ L/ max(u t h). However, this time is 
expected to be very small since max(u t h) is likely to much larger than the sound speed. 
Therefore, the reduction of exactly pollutant-free fraction, P(t), by molecular diffusion 
may be considered as being essentially instantaneous. 

For our astrophysical applications, we need the fraction, P(Z c ,t), of the flow with Z 
below a finite critical value, Z c , rather than the exactly pristine fraction. Obviously, it 
takes finite time for molecular diffusivity to enrich all the fluid elements in the system 
up to a finite threshold, Z c . In fact, during a short time interval, the degree of pollution 
by molecular diffusivity alone is negligible even at small distances from the pollutant 
source, and the entire system is practically unpolluted. Therefore, the observation of the 
rapid/immediate erasure of exactly pristine gas by molecular diffusivity is not directly 
relevant to the astrophysical problem of primordial star formation. 

Because k is usually tiny in practical environments, such as in the interstellar media 
of galaxies, enriching all the fluid elements to a concentration level of, say, Jb 10~ 8 , 
by the molecular diffusivity alone is very slow (see discussions in Pan and Scalo 2007). 
The presence of a turbulent velocity field greatly accelerates the mixing process, making 
the reduction of P(Z c ,t) much faster. We find that the timescale for the reduction of 
P(Z c ,t) with a small Z c is basically determined by the rate at which the turbulent 
stretching produces small-scale structures and is essentially independent of k. 

4.1. The Mapping Closure Model 

We calculate the fraction P(Z c ,t) predicted by the mapping closure model. From eq. 
(|3.7[) . it is straightforward to find that 




(4.2) 



— OO 
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where the upper limit rj c (t) satisfies X(r) c (t);t) — Z c . For a given value of Z c , the limit 
rj c (t) changes with time as the mapping function evolves, and for our initial bimodal PDF 
with two delta functions, rj c (t) can be computed using eq. (|3.8[) . 

From that equation, we see that Z c = corresponds to r/ c — oo at all times after 
t = 0. Therefore, the mapping closure model predicts that Pit) is zero at any time 
t > 0, or that the fraction of exactly pollutant-free mass deceases to zero instantaneously. 
This is consistent with our discussion above that the molecular diffusivity alone tends 
to immediately remove fluid elements with exactly zero concentration. The mapping 
closure model inherits this particular property of molecular diffusion, because the effect 
of diffusivity as a Laplacian term is treated directly. The model destroys the initial delta 
function at Z = instantaneously. However, this does not suggest that p(0; t) becomes 
finite immediately. At the early evolution stage, p(Z\ t) does have an infinite peak at 
Z = 0, but the peak is less singular than a delta function (S(Z)) in the sense that 
J"<f p(Z';t)dZ' -> in the limit Z ->• 0. 



4.2. The Nonlinear Integral Models 

Unlike the mapping closure model, the nonlinear integral models preserve the singularities 
at Z = and Z = 1. More specifically, the amplitudes of the delta functions at Z = and 
Z = 1 decrease with time, but they are never completely destroyed, such that exactly 
pollutant-free mass can survive in these models, and P(t) remains finite at any finite time. 
This is inconsistent with our earlier observation that the molecular diffusivity tends to 
reduces P(t) to zero immediately. The reason is that the effect of molecular diffusivity 
is not incorporated directly in these models, instead it is included implicitly through 
the function J(Z; Zx, Z 2 ). Despite the inconsistency, we find that the integral models 
are useful for understanding the pollution of fluid elements with very low (but nonzero) 
concentration by turbulent mixing. Below we derive an equation for the fraction, P(t), 
of exactly pollutant-free mass from these models. 

We consider the general model represented by eq. (|3.10|) . Integrating this equation in 
the range [0, Z] and taking the limit Z — > 0, we have, 



dP(t) 
dt 



7(t) J dZ lP (Z i; t) J dZ 2 p(Z 2 ;t) lim J dZ'j(Z';Zi,Z 2 )-P(t) . (4.3) 



The last integral in the triple-integral term in the limit Z — > can be written as 

f J(Z'; Zi, Z 2 )dZ' where + represents the upper integral limit approaching zero from 
the positive vicinity. We first note that this integral is zero if both Z\ and Z 2 are pos- 
itive because J(0; Z\, Z 2 ) — for Z\ > and Z 2 > (see §3.3). We next assume that 
J(Z; Z\, Z 2 ) at Z = Z\ and Z — Z 2 is nonsingular or less singular than a delta function 

for Z\ 7^ Z 2 (meaning that jj^ 1 J(Z; Z%, Z 2 )dZ = and f z - dJ(Z; Z\, Z 2 )dZ = 0, where 

Z\ < Z 2 is assumed without loss of generality). This assumption is clearly satisfied for 
Curl's model and the model with uniform J(Z; Zi, Z 2 ) (eq. (| 3 . 1 lj) V With this assump- 
tion, it is straightforward to see that J Q J{Z'; Z\, Z 2 )dZ' is finite only if both Z\ = 

and Z 2 = 0. In that case, we have J Q ° dZ' J{Z'\ 0, 0) = 1 from the normalization of 
J(Z; Zi, Z 2 ). This observation suggests that the contribution to the triple integral in eq. 
(|4.3[) comes only from Z\ and Z 2 values in an infinitesimal range around zero. With such 
infinitesimal ranges of Z\ and Z 2 , the first two of the three integrals contribute factors 

of J Q p(Z\\t)dZ\ and L p(Z 2 ;t)dZ 2 , respectively. As both these factors are equal to 
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P(t), we have the following equation for P(t), 

d_m__ m-p) , 44) 

dt T int ' 1 ' 

where r; n t = 7(t) _1 is used for the convenience of notations. If the mixing timescale T m 
is constant, the variance decay requirement gives r in t = T m /2 or r m /3 for Curl's model 
and the model with uniform J(Z; Z\, Z2), respectively (see §3.3). 

The equation gives an interesting physical picture for mixing of the unpolluted mass 
in turbulent flows: the pristine fraction is reduced when turbulent stretching brings the 
pollutant-free fluid elements (with a fraction of P(t)), and the rest of the flow (with a 
fraction of 1 — P(t)), which has been polluted by sources or previous mixing events, close 
enough for molecular diffusivity to homogenize (Pan and Scalo 2007). 

Eq. (|4.4p has a simple analytic solution, 

p (t) = r^r> ( 4 - 5 ) 



P + (1-P )exp(^) 



where Pq is the initial fraction of unpolluted mass, and we have assumed Ti n t is constant 
with time. Although it is derived for the fraction of exactly pollutant-free mass, we will 
show in §5 that, in certain physical regimes, this equation can be used to fit our numerical 
results for P{Z c ,t) with a finite threshold Z c . 

4.3. The Self- convolution Models 

The self-convolution models introduced in §3.4 also preserve the initial singularities at 
Z = and Z = 1, since they are essentially extensions of Curl's model. Again we derive 
the equations for the fraction, P(t), of exactly pollutant-free mass from the convolution 
models, which will be used later to understand the mass fraction of nearly-pristine, but 
Z / 0, flow. We first decompose the concentration PDF as 

p(Z;t)=P(t)S(Z)+p c (Z;t), (4.6) 

where p e (Z; t) denotes the concentration PDF in the enriched/polluted part of the flow, 
and it satisfies that lim^^o J Pc(Z'; t)dZ' = 0. The Laplace transform of eq. (|4.6p gives, 

p((;t)=P(t)+p e (C,t), (4.7) 

where p c (£;i) is the Laplace transform of p e (Z;t). In the limit £ — > +00, p c (£;t) ap- 
proaches zero because limz->.o J Pe{Z'\t)dZ' = 0. 

Inserting eq. (|4.7p to eq. (|3.13|) for the model of Venaille and Sommeria (2007), and 
taking the limit ( — > +00, we find, 



dP{t) _ Pln(P) 

dt 7"con 



(4.8) 



where r con = 7 1 , and we used the fact that p e ((;t) — > and (d<;p c ((;t) ~J> as ( 
approaches infinity. If r con is constant, the equation is solved by, 

P(t) = P exp(t/T - ) , (4.9) 

which can also be obtained from the solution for p(£; t) given in Venaille and Sommeria 
(2007). We will show that eq. (|4.9p provides a useful fitting function for our simulation 
data for P(Z c ,t) with finite Z c in a transonic flow. 

Similarly, we can derive an equation for the pristine fraction from the generalized 
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version, eq. ()3.14p . of the self-convolution models, 

dt t 

Assuming both n and r con are constant with time, the solution of the equation is, 

Po 



P(t) 



P, 



l/n 



(1 - P 1/n ) exp (i/r coa ) 



(4.10) 



(4.11) 



For n = 1, the equation reduces to eq. (|4.5|) for the nonlinear integral models, and in the 
limit of n — > oo, it approaches eq. (|4.9p for the continuous convolution model. We will 
use eq. (|4.1ip to fit our simulation results for scalars in a highly supersonic flow, taking 
T con and n as fitting parameters. 



5. Numerical Simulations 

To test the theoretical models and fix their parameters, we carried out numerical simu- 
lations for mixing in hydrodynamic turbulent flows using the FLASH code (version 3.2), a 
multidimensional hydrodynamic code (Fryxell et al. 2000) that solves the Riemann prob- 
lem on a Cartesian grid using a directionally-split Piecewise-Parabolic Method (PPM) 
solver (Colella & Woodward 1984; Colella & Glaz 1985; Fryxell, Midler, & Arnett 1989). 
We evolved the hydrodynamic equations, 

§ + V-(pv)=0, 

| + v.W = ^|f, (5.1) 
ot p 

on a 512 3 grid for a domain of unit size with periodic boundary conditions. We adopted an 
isothermal equation of state, p = pC s , with a constant sound speed, C s . The isothermal 
equation of state is commonly used to imitate the nearly constant temperature in some 
interstellar environments, and is a convenient assumption to investigate the effects of 
compressibility in interstellar turbulence. Our code does not explicitly incorporate a 
viscosity term, and the kinetic energy is dissipated by numerical diffusion. A large-scale 
solenoidal external force, f , was applied to drive and maintain the turbulent flows. This 
driving force was taken to be a Gaussian stochastic vector with an exponential temporal 
correlation function. We generated f in Fourier space and included all independent modes 
with wave numbers in the range from 2ir and 67r. Each independent mode was given the 
same amount of power. We defined a forcing length scale as Lf = J ^-Pf(k)dk/ J Pf(k)dk, 
with Pf(k) being the power spectrum of the driving force, and found that Lf was equal 
to 0.46 box size for our driving scheme. 

We adjusted the amplitude of the driving force to obtain a transonic flow with rms 
Mach number M = 0.9 and a supersonic flow with M = 6.2. We refer to the two flows 
as flow A and flow B. The rms Mach number was defined as the density- weighted rms 
velocity, v lms , divided by the sound speed, C s , and was computed from the temporal 
average after the flow reached a statistical steady state. We defined a flow dynamical 
timescale as Td yn = L{/v TmB . The simulation setup for the turbulent flows is the same as 
that in Pan & Scannapieco (2010), to which we refer the interested reader for details. 

To study mixing, we solved the advection equation for a number of decaying scalars, 
which were added to the flow once the turbulence had become fully developed and sta- 
tistically stationary. The initial concentration field of the decaying scalars was bimodal, 
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consisting of pure pollutants and completely unpolluted flow. The initial pollutant region 
was chosen to be a single cube located right at the center of the simulation box. Within 
this cube, we set the concentration field, C to be unity, i.e., the flow material there was 
taken be pure pollutants, and outside of the cube we set C = 0, i.e., the flow there 
was completely pollutant free. This initial condition was chosen for its simplicity, and 
it suffices for the purpose of illustrating the general problem and testing the theoretical 
models. 

An important parameter for the initial condition is the pollutant fraction, Pi, i.e., the 
ratio of the pollutant mass to the total mass in the simulation box. Clearly, the fraction 
fixes the initial pristine fraction, Pq = 1 — Pi. We considered three scalars in each of 
our flows and set the initial pollutant fraction to be Pi = 0.5, 0.1 and 0.01, respectively. 
In the M = 0.9 flow, we name the three scalars with Pi— 0.5, 0.1, and 0.01 as Al, 
A2 and A3, respectively. The corresponding cases in the M — 6.2 flow are named Bl, 
B2 and B3. The exact values for Pi were achieved by tuning the size of the pollutant 
regions. Smaller values of Pi would be also of interest for mixing of heavy elements in 
the interstellar media of early galaxies. However, for Pi <C 0.01, the size of the pollutant 
region becomes smaller than the integral scale of our simulated flows. This gives rise to 
complications in the evolution of the unpolluted (or slightly polluted) fraction. Smaller 
initial pollutant fractions will be investigated in a followup study. 

Similar to the case of kinetic energy dissipation, the scalar dissipation (or homogeniza- 
tion) is also through numerical diffusion in our simulations. The diffusion scale is thus 
close to the resolution scale. To examine whether our results depend on the amplitude 
of numerical diffusion, we performed the same runs at a lower resolution, 256 3 , and con- 
ducted a convergence study. We found that the timescale for the evolution of P(Z c ,t) 
with Z c ~ 10~ 8 already converged at the resolution 512 3 . 



6. Results 

6.1. The Concentration Field 

In Fig. [TJ we plot the evolution of the concentration field of scalar A2 in our simulated 
flow with M = 0.9. The four panels correspond to the log of the concentration field on 
a slice (z — 0.5) of the simulation grid at four snapshots with t =0.12, 0.5, 0.9 and 1.5 
Tdyn, respectively. The color table is in logarithmic scale and the lower limit was chosen 
to be 10 -8 , so that the part of the flow with concentration below a small threshold, Z c , 
is visible in the figure. The size of initial pollutant at the center of the box was set to 
be 0.47 in units of the box size, such that the initial pollutant fraction is 0.1. With time, 
the turbulent flow transports and spreads out the pollutants, exposing them to more and 
more pristine fluid elements. Turbulent stretching by vortices and shear continuously 
produces concentration structures at smaller and smaller scales. Pristine fluid elements 
are contaminated when encountering a pollutant/polluted structure within a distance 
smaller than the diffusion scale. At t = 1.5rd yn , almost the entire flow is polluted, and 
the mass fraction of the flow with Z ^ 10 -8 becomes negligibly small. Cliff structures 
typical of passive scalar fields advected in incompressible turbulence are clearly observed 
in panels (c) and (d). 

Fig.[Ushows the concentration field of scalar B2 in our M — 6.2 flow. At the four snap- 
shots selected here, the density- weighted concentration variances are close to those for 
scalar A2 at the corresponding snapshots in Fig. Q] It appears that, at similar variances, 
the unpolluted volume in the M = 6.2 flow is significantly larger than in the M = 0.9 
case. The existence of strong compressions and expansions in a highly supersonic flow 
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(a) (b) 




Figure 1. Log of the concentration field of scalar A2 on a slice of the simulation grid at snapshots 
with t = 0.12 (a), 0.5 (b), 0.9 (c) and 1.5 T dyn (d). The scalar is advected in the M = 0.9 flow. 
The size of the initial pollutant cube is 0.47 box size. The color table ranges from 10 -8 to 1, 
with the white color representing regions with concentration Z 10 -8 . 



gives rise to a very different geometry for the scalar field. Most of the volume in a highly 
compressible flow is occupied by expanding events, and the flow expansion tends to pro- 
duce more coherent scalar structures, as a passive scalar follows the expansion. Therefore, 
scalar B2 appears to be smoother than A2 in Fig. [TJ The edge-like structures in Fig. [5] 
are produced by the compression of velocity shocks, which amplifies the scalar gradient 
across the shocks. Although the visual impression of Fig. [2] is dominated by the effect 
of compressible modes, it is actually the solenoidal modes that contain the majority of 
kinetic energy in the flow and provide the primary contribution to the scalar cascade 
even at high Mach numbers (Pan & Scannapieco 2010). The interested reader is referred 
to Pan and Scannapieco (2010, 2011) for detailed discussions of scalar structures as a 
function of the flow Mach number and the relative role of solenoidal and potential modes 
for mixing in supersonic turbulence. 



6.2. The PDF Evolution 

Fig. [3] plots the PDFs of scalars A2 (left panel) and B2 (right panel) as a function of 
time. The two scalars were evolved in our simulated flows with M = 0.9 and M = 6.2, 
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Figure 2. Same as Fig. [TJ but for scalar B2 in the M = 6.2 flow. The four snapshots, (a), (b), 
(c) and (d), correspond to t = 0.11, 0.65, 1.1 and 1.7 Td yn , respectively. The size of the initial 
pollutant cube is also 0.47 box size, and the initial pollutant fraction is 0.1. 

respectively. The initial pollutant fraction, Pi, is 0.1 for the two scalars, meaning that the 
amplitude of the initial spike at Z = is 9 times higher than that at Z = 1 . Turbulent 
mixing reduces the heights of the two spikes, and gradually fills the concentration space 
in between. Eventually a central peak forms around the average concentration, and then 
the PDF narrows down toward the peak. 

The lines in Fig. [3] are the prediction of the mapping closure model. At each time, the 
predicted PDF has the same value of variance as that from the simulation (data points). 
This is equivalent to properly choosing the timescale \^{t)/n in eq. (|3.6[) so that the 
variance evolution from the model matches the simulation result. The model prediction 
is in good agreement with the data at the central part and the right (high-Z) tail of 
the scalar PDF. However, the mapping closure considerably underestimates the left PDF 
tail at intermediate to late times. A detailed discussion of the discrepancy between the 
prediction of the Gaussian mapping closure and simulation results at late times is given 
in Girimaji (1992). The weakness of the mapping closure is also discussed by Duplat and 
Villermaux (2008). It appears that, for scalar A2, the agreement between the mapping 
closure and the simulation data becomes better in the long time limit. 

In the right panel of Fig. [3] for scalar B2 in the M — 6.2 flow, we see that at early 




Figure 3. PDF evolution of scalars A2 (a) and B2 (b) in the M = 0.9 and M = 6.2 flows, 
respectively. The initial pollutant fraction, Pi, for the two scalars is 0.1. Lines are the predicted 
PDFs from the mapping closure model with the same values of variance as in the simulations. 



times there is also an acceptable agreement between the mapping closure prediction and 
the simulation results. At later times, the discrepancy at the left tails between the model 
and the data is larger than the M = 0.9 case. This is because the PDF tails broaden with 
increasing Mach number, as previously found in the simulations of Pan & Scannapieco 
(2010). The origin of this effect was argued to be related to the increasing degree of 
intermittency of the velocity field as the Mach number increases. 

The mass fraction of unpolluted or slightly polluted flow with Z <S 10~ 8 — 10~ 5 cor- 
responds to the far left tail of the PDF with Z well below the minimum value shown 
in Fig. [3J Therefore, Fig. [3] does not contain direct information for the part of the PDF 
of our primary interest. Nevertheless, the left PDF tails shown in Fig. [3] imply that the 
mapping closure model is likely to significantly underestimate the unpolluted fraction at 
intermediate to late times especially for scalar B2. The expectation is confirmed in §6.4. 

In the left panel of Fig. we compare the prediction of the nonlinear integral model 
with uniform J(Z; Z±, Zi) to the simulation data for scalar A2 in the M = 0.9 flow. 
The performance of this model for the PDF evolution is poor. At very early times, the 
predicted PDF appears to be flat in between the initial spikes, reflecting a "memory" 
of the uniform function J(Z; Z±, Z\). As mentioned earlier, a problem of the nonlinear 
integral models is that they predict excessively fat tails in the long time limit. This 
is seen in the left panel of Fig. 21 which shows that at large t the model significantly 
overestimates the left tails. We find that at late times the nonlinear integral model also 
overestimates the PDF tails for scalars in the M — 6.2 flow (not shown). Although the 
nonlinear integral models do not give good predictions for the PDF evolution, we find 
that they provide useful fits to the pristine fraction in certain physical regimes (see §6.4). 

In the right panel of Fig. gl we fit the PDF of scalar A2 in the M = 0.9 flow with 
the Gamma distribution, eq. (|3.16l) . predicted by the model of Villermaux and Duplat 
(2003). For each line, the value of n is chosen such that the variance of the Gamma 
distribution is equal to that from the simulation. The scalar PDF at t ^ lTd ym does not 
have a Gamma distribution shape, and is not shown in this panel. At the four selected 
times in the figure, however, the Gamma distributions fit the simulation data quite well. 
The agreement is significantly better than the mapping closure for t between 1.8 and 2.6 
Tdyn (corresponding to 1.1 ^ n ^ 4). We find that n(t) increases exponentially with time, 




Figure 4. PDF evolution of scalar A2 in the M — 0.9 flow. Lines are predictions of two 
models, (a): the nonlinear integral model with uniform J(Z; Z\, Z2) (eq. p. lip ), (b): Gamma 
distributions as predicted Villermaux and Duplat (2003). 



corresponding to the exponential decay of the scalar variance at late times (see §6.3), as 
the variance of the Gamma distribution, eq. (|3.16p . goes like (Z) 2 jn. This is in contrast 
to the experimental result, n(t) oc t 5 / 2 , found by Villermaux and Duplat (2003). The 
reason for the difference is that our simulated flows are maintained at a steady state by 
a driving force and are statistically homogeneous, while the experiments by Villermaux 
and Duplat (2003) are for decaying flows dominated by a mean shear. An exponential 
decay is also found in a sustained flow by Villermaux et al. (2008). 

We point out that the continuous convolution model of Venaille and Sommeria (2007) 
predicts that, if the scalar PDF is given by a Gamma distribution at a given time, then 
the PDF will remain a Gamma distribution at all subsequent times. Therefore, if one 
starts to use the continuous convolution model at a time when the scalar PDF has evolved 
to a Gamma distribution, its prediction for later times would be the same as the model 
of Villermaux and Duplat (2003). However, unlike Villermaux and Duplat (2003), the 
model of Venaille and Sommeria (2007) does not predict that the Gamma distribution 
is an attractive solution that the scalar PDF always reaches at the late evolution stage. 

For scalar A2 in the M = 0.9 flow, the model of Villermaux and Duplat (2003) starts 
to apply at t = 1 — 2r ( i yn . At this time, the pristine mass fraction has already decreased 
to very small values, and thus the model is not suitable to study the pristine fraction. 
We also tried to fit Gamma distributions to the scalar PDFs in the M = 6.2 flow, and 
found they underestimate the left tails, which are broader than the M = 0.9 case. 

In summary, we found that in the M = 0.9 flow the mapping closure gives acceptable 
fits to the scalar PDF at early times, but significantly underestimates the left PDF tails at 
late times. Starting from ~ 1.8rd yn , the scalar PDF is better fit by a Gamma distribution, 
which is predicted by the model of Villermaux and Duplat (2003). Since all the models 
we considered were originally developed to study mixing in incompressible turbulence, 
they were not expected to perform well in highly compressible flows. In the M — 6.2 flow, 
the left PDF tail is broader, and no models were found to give satisfactory predictions 
for the scalar PDF at late times. 

6.3. The Variance Decay 

In Fig. [5l we show the evolution of the density- weighted concentration variance for scalars 
A2 and B2 in the M = 0.9 and M = 6.2 flows, respectively. In the left panel, the time is 




Figure 5. Evolution of the density-weighted concentration variance for scalars A2 and B2 in 
M = 0.9 and M = 6.2 flows, respectively. The variance is normalized to its initial value, (a): 
time normalized to the flow dynamical timescale. (b): time normalized to the acoustic timescale. 



normalized to the flow dynamical time Td yn , while in the right panel it is normalized to 
the acoustic time r ac defined as Lf/C s , the time for the sound wave to cross the driving 
scale of the flow, Lf. From the definition, we have r ac = MTd yn . Therefore, the variance 
decay curves shift to the right by a factor of 1.1 for scalars A2 and to the left by a factor 
of 6.2 for B2, when the normalization timescale is switched from Td yn to r ac . From Fig. 
[5j it is clear that the flow dynamical time is a more relevant timescale for the scalar 
variance decay in supersonic isothermal turbulent flows. 

As seen in the left panel of Fig. [5j at early times the variance decreases slowly, cor- 
responding to the initial development of scalar structures toward small scales. In this 
transient period, the parameter "/(t) in the models presented in §3.3 and §3.4 would be 
time-dependent if the model is required to match the scalar variance decay. At t ^ 0.5Td yn , 
the variance decays exponentially, and "/(t) would be constant. When normalized to the 
dynamical timescale, the variance decay of scalar A2 is slightly faster than B2, and the 
decay timescale is measured to be r m = 0.61 and 0.73rd yn , for A2 and B2, respectively. 
These results are consistent with Pan & Scannapieco (2010), who found that the variance 
decay timescale in units of the flow dynamical timescale has a weak dependence on the 
flow Mach number, increasing by <^ 20% as M goes from 1 to 6. This slight increase is 
due to the fact that compressible modes are less efficient at enhancing the mixing rate 
than solenoidal modes (Pan & Scannapieco 2010). 

We also measured r m for the other four scalars included in our simulations, and found 
that, in each flow, the mixing timescales of the three scalars are close to each other. The 
timescale for the third scalar (i.e., A3 or B3) is slightly smaller than the other two. On 
average, the mixing timescale is r m ~ 0.6rd yn for the three scalars in the M = 0.9 flow. 
In the M = 6.2 flow, the average mixing timescale is ~ 0.7rd yn - 

6.4. The Fraction P(Z c ,t) 

Next we computed the mass fraction, P(Z c ,t), of fluid elements with Z smaller than 
different thresholds, Z c , from the simulation data. We found that the flow with exactly 
zero concentration (i.e., the special case with Z c = 0) is erased rapidly by numerical 
diffusion. This is because in each time step the unpolluted computation cells adjacent 
to those with nonzero Z obtain a finite (but tiny) concentration value due to numerical 




Figure 6. Mass fraction of fluid elements with Z ^ 10 for scalars Al and A2 in the M = 0.9 
flow. Dashed lines correspond to the prediction of the mapping closure model. Solids lines are 
the best fits using the continuous convolution model with r con = 0.35 and 0.37 Td yn for Al and 
A2, respectively. The dotted lines are fits by eq. (14. 5|) from the nonlinear integral models. 



diffusion. Therefore, after a small number of time steps, no exactly pollutant-free cells 
were left in the simulation box. However, in such a short time, the degree of pollution in 
most cells due to numerical diffusion is extremely small, and the concentration level in 
these cells was much smaller than any threshold of practical interest. The rapid pollution 
of completely pristine mass by numerical diffusion in the simulations is similar to the effect 
of molecular diffusivity, which tends to reduce P(t) to zero instantaneously, although the 
numerical diffusion probably has a different form and a much larger amplitude than the 
realistic molecular diffusivity. 

The threshold of interest for astrophysical applications is Z c <; 10 -8 . For this finite 
threshold, the timescale at which P(Z c ,t) decreases is significant, and is on the order of 
the flow dynamical time. A comparison of the same simulation runs at two resolutions, 
256 3 and 512 3 , shows that the timescale for P(Z C , t) with Z c ~ 10~ 8 is independent of the 
amplitude of numerical diffusion. These suggest that the reduction rate of P(Z c ,t) with 
Z c 10~ 8 is mainly determined by the large-scale properties of the flow. The behavior 
of P(Z c ,t) as a function of time is similar for different values of Z Cl given that Z c is 
much smaller than the average concentration. The evolution timescale of P(Z c ,t) only 
has a weak logarithmic dependence on Z c . Below we only present results for Z c = 10~ 8 . 
Similar results are found for Z c in the range from 10~ 8 to 10 -5 . 

6.4.1. The M = 0.9 Flow 

In Fig. [SI we plot the mass fraction of fluid elements with Z ^ 10 -8 as a function of time 
for scalars Al and A2 in the M = 0.9 flow. The initial pollutant fraction, Pi, is 0.5 and 
0.1 for the two scalars, respectively. The data points are results from the simulations. The 
dashed lines correspond to the prediction of the mapping closure model. For this model, 
the fraction, P(Z C , t), at a given time is calculated from the predicted PDF with the same 
scalar variance as that in the simulation. In Fig. [6j we see the mapping closure model is 
in good agreement with the data points for scalar Al. However, the model prediction is 




Figure 7. Mass fraction of fluid elements with Z < 10 for scalar A3 in the M = 0.9 flow. 
The dashed line corresponds to the nonlinear integral model with Tint = 0.24rd yn . The solid line 
is the best fit obtained by combining the nonlinear integral model for the early phase and the 
continuous convolution model for the later phase. The timescales used in the two phases are 
Tint = 0.24Td yn and r con = 0.36rd yn , respectiely. The inset shows the same data points and lines, 
but with the vertical coordinate on a linear scale. 



well below the data points for scalar A2. This was expected from the left panel of Fig.[3l 
which shows that the mapping closure model underestimates the left PDF tail of scalar 
A2 at intermediate to late times. We found that the model also underestimates P(Z c ,t) 
for scalar A3, which had an initial pollutant fraction of 0.01, and the discrepancy is even 
larger than the case of A2. 

The solid lines in Fig.[5]are from the continuous convolution model of Venaille and Som- 
meria (2007), i.e., eq. (|4.9[) . and they match the data points quite well. The timescales r con 
used in the fits are 0.35 and 0.37 Td ym , respectively, for scalars Al and A2. As discussed 
in §4.3, eq. (|4.9|) was originally derived for the fraction, P(t), of exactly pollutant-free 
mass. The good agreement between eq. (|4.9|) and the simulation data shows that the 
model actually provides an excellent fitting function for the fraction, P(Z c ,t), with a 
finite (but small) threshold Z c . It also suggests that the continuous convolution process 
is a good physical description for the erasure of unpolluted (or slightly polluted) flow by 
turbulent mixing, if the initial pollutant fraction, Pi, is larger than ~ 0.1. 

The dotted lines in Fig. [6] shows the fits using eq. (|4.5[) from the nonlinear integral 
model. The timescale T; n t was chosen to be 0.28rH yn and 0.3-7d yn for scalars Al and A2, 
respectively. This model predicts an exponential decrease at late times, which is much 
slower than found in the simulation. The overestimate is because the model produces 
excessively broad PDF tails in the long time limit (see left panel of Fig. 2]). 

Fig.[7]shows the evolution of P(10 _8 ,t) for scalar A3, whose initial pollutant fraction 
Pi = 0.01. The inset plots the same data points and model fits, but the y-axis is on a 
linear scale. Unlike the case of scalars Al and A2, the data points for scalar A3 cannot 
be well fit by the continuous convolution model with a single timescale right from the 
beginning. In fact, P(10~ 8 ,t) exhibits different behaviors at early and late times. 

The early phase can be well fit by eq. (|4.5jl from the nonlinear integral model with 
r; nt = 0.24rd yn . This is shown as a dashed line in Fig. and from the inset we see 
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the line matches the simulation data over an extended time range before P(lCP 8 ,i) 
decreases to ~ 0.3. The timescale Tint used here corresponds to ~ 0.4r m since the variance 
decay timescale r m for scalar A3 is about 0.6Td yn . This is in between 1/2 and 1/3 r m , 
the expected values of Ti n t for Curl's model and the model with uniform 3(Z\ Z\, Z%), 
respectively (see §4.2). The dashed line starts to significantly overestimate the simulation 
results when P(10 -8 , t) becomes smaller than ~ 0.3. Again, this is because the nonlinear 
integral models significantly overpredict the PDF tails at late times. 

We find that the late-time behavior of P(10 -8 , t) can be well described by the contin- 
uous convolution model. In fact, this model starts to give a satisfactory fit quite early, 
right after P(10~ 8 , t) becomes smaller than ~ 0.8. The best-fit value of the timescale T ccm 
is 0.36rd yn , which is almost the same as the values used to fit the data for scalars Al and 
A2. This, together with the results for scalars Al and A2, suggests that the continuous 
convolution model applies if the mass fraction of pollutants or the polluted flow is larger 
than 0.1 — 0.2. We point out that there is an extended range of P(10~ 8 ,i) (from 0.3 to 
0.8), where both the nonlinear integral model and the continuous convolution model can 
well match the data. 

We give a physical speculation for why the early phase of scalar A3 is better fit by 
the nonlinear integral model than the continuous convolution model. For this scalar, the 
amount of pollutants or polluted mass is small at early times. The limited availability of 
pollution sources leads to a relatively low frequency of contact between the polluted and 
unpolluted flow. As a consequence, the pollution process would involve less convolutions 
at this stage, and thus may be better captured by the nonlinear integral model, which can 
be roughly viewed as a discrete self-convolution model in Laplace space. The simulation 
results presented above suggest that the mixing events between the unpolluted and the 
polluted fluid elements become frequent enough to trigger the continuous convolution 
process, when the polluted fraction is larger than 0.1-0.2. 

The solid line in Fig. [7] is obtained by combining the best fits for the early and late 
phases using the nonlinear integral model and the continuous convolution model, re- 
spectively. Clearly, this line is in excellent agreement with the simulation results. We 
connected the two phases at time, £0.5, when P(10 _8 ,i) = 0.5, i.e., the second phase is 
fit by 0.5 cxp ^*~*° ' 5 '/ Tcon L As mentioned earlier, the best-fit timescales for the two phases 
are T ln t = 0.24rd yn and r con = 0.36rd yn , respectively. The choice of connecting the two 
phases at P(10~ 8 ,t) = 0.5 is somewhat arbitrary. In fact, combining the two models at 
any time with P(10 _8 ,i) between ~ 0.8 and ~ 0.3 yields satisfactory results. 

The timescales, r; n t and r con , were set to be constant in all our fits to the simulation 
results for P(Z C , t). These timescales can be a function of time in general. In fact, from the 
consideration of the scalar variance decay, these timescales would be time-dependent at 
early times when the variance decay is slower than exponential, and then become constant 
at t J> 0.5rd yn when the exponential decay starts (see §6.3). The reason why assuming 
constant timescales applies perfectly for the evolution of the pristine mass fraction, but 
not for the scalar variance decay at all times is probably that the pollution of pristine 
flow is physically simpler. A fast variance decay relies on the full development of scalar 
structures toward the diffusion scale, and that is why the decay is slow at early times as 
the cascade just starts. On the other hand, the pollution of pristine mass does not need 
to wait for the scalar structures to be fully developed at small scales: the pollution occurs 
whenever the unpolluted fluid elements are brought into contact with the pollutants or 
the polluted flow. This happens right away once the pollutants are released into the flow. 
This is perhaps why using constant timescales, Tint and r con , right from the beginning 
could give successful fits to the evolution of P(Z c ,t). 
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Figure 8. Mass fraction of fluid elements with Z < 10 for scalars Bl, B2 and B3 in the 
M = 6.2 flow. The solid lines for Bl and B2 correspond to the fits by the continuous convolution 
model (eq. I4.9|l with r con = 0.46rd yn for both cases. The line for B3 combines the nonlinear 
integral model with Tint = 0.30rd yn for the early phase and the continuous convolution model 
with T con = 0.5lTd yn for the later phase. 



6.4.2. The M = 6.2 Flow 

Fig.|8]shows our results for scalars Bl, B2 and B3 with initial pollution fractions of 0.5, 
0.1, and 0.01, respectively, in the M = 6.2 flow. We find that the mapping closure model 
significantly underestimates P(10 _8 ,i) for all the three scalars, and the discrepancy is 
much larger than the case of the M = 0.9 flow. In Fig.0 we attempt to fit the simulation 
data for scalars Bl and B2 with the continuous convolution model [eq. (|4.9|) ] of Venaille 
and Sommeria (2007), as in the M = 0.9 case. The timescale, r con , used in the fitting 
lines is 0.46 r^yn for both Bl and B2. Again the fitting curve for scalar B3 consists of two 
phases that connect at P(Z c ,t) = 0.5. The early phase is fit by the nonlinear integral 
model with r; n t = 0.3THy n , and the later phase uses the continuous convolution model 
with r con = 0.51rd yn . The parameter choice here gives priority to satisfactorily matching 
the data points at early times. 

All the timescales chosen in Fig. [5] are larger than the corresponding values used for 
scalars in the M = 0.9 flow. This is caused by two effects. First, as mentioned earlier, 
when normalized to the flow dynamical time, the variance decay timescale in the M = 6.2 
flow is slightly larger than in the M = 0.9 case. Second, as shown in Fig. [3J the left tail 
of the scalar PDF broadens with increasing Mach number of the advecting flow. This 
means that, with the same concentration variance, the scalar PDF in the M = 6.2 flow 
contains a larger probability at low concentration levels. Both effects tend to result in a 
slower decrease of P(Z c ,t) in the flow with higher M. The second effect appears to be 
stronger than the first one, and it also explains why the same models that well match 
the results in the M = 0.9 flow significantly underestimate P(Z C , t) in the M = 6.2 flow 
at late times (see Fig. [8]). The fitting quality of the lines in Fig. [8] is good at early times 
when P(Z c ,t) J> 0.1, and generally acceptable before P(Z c ,t) decreases to 0.01. Below 
that, however, significant discrepancy appears. 

We find that the simulation results for Bl and B2 in the M = 6.2 flow can be better fit 
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Figure 9. Mass fraction of fluid elements with Z < 10 for scalars Bl, B2 and B3 in the 
M = 6.2 flow. The data points are the same as in Fig. [8] Eq. (|4.11fl with n — 4.6 is used to 
match the simulation data. The timescale, r con , in the equation is set to 0.40 and 0.41rd yn for 
scalars Bl and B2, respectively. The line for B3 is a combination of the nonlinear integral model 
with Ti n t = 0.30rd yn for the early phase and eq. (|4.1ip with r con = 0.42rd yn and n — 4.6 for the 
later phase. 

by eq. (|4.11|) from the generalized self-convolution model (see eq. (|3.14l) in §3.4; Duplat 
& Villcrmaux 2008). The parameter n in this equation controls the shape of the fitting 
curve. Fig. |9] shows our results using this equation to fit the simulation data. The data 
points here are the same as in Fig. [8] Eq. (|4.11[) with n = 4.6 can well fit the simulation 
data for scalars Bl and B2 at all times and for scalar B3 in the late phase. For Bl and 
B2, the best-fit timescale r con is, respectively, 0.40 and 0.41 Td yn . For the early phase of 
scalar B3, we used the same nonlinear integral model as in Fig. [8] with Ti nt = 0.3rd yn . 
The late phase is fit by eq. (|4.11j) with T con = 0.42rd yn and n = 4.6. We combined the two 
phases at P{Z c ,t) = 0.5. A comparison of Fig. |8] and Fig. [9] shows that, with eq. (|4.11|) 
from the generalized convolution model, the fitting quality is significantly improved. 

We point out that eq. (|4.11l) is used simply as a fitting function. The generalized 
convolution model (§3.4) behind this equation does not address the effects of shocks and 
compressibility on the PDF of passive scalars in supersonic turbulence. There is thus 
no physical reason why it provides successful fits to the pristine mass fraction in the 
M = 6.2 flow. A physical closure model is motivated to successfully explain and match 
the evolution of P(Z c ,t) in highly supersonic turbulence. 

7. Conclusions 

Motivated by the process of primordial star formation in the first generation of galaxies, 
we investigated the general problem of how the unpolluted flow material in compressible 
turbulence is contaminated by mixing. We approached this problem using both theo- 
retical modeling and numerical simulations. The theoretical approach is based on the 
probability distribution method for turbulent mixing, since the fraction of the unpol- 
luted or slightly polluted mass corresponds to the left tail of the concentration PDF. We 
first derived an equation for the concentration PDF with density- weighting, where the 



2G 



L. Pan, E. Scannapieco & J. Scalo 



advection term exactly conserves the global PDF. We then considered several existing 
closure models for the diffusivity term in the PDF equation, including the mapping clo- 
sure model (Chen et al. 1989), the nonlinear integral models (Curl 1963, Dopazo 1979, 
Janicka et al. 1979) and the self-convolution models (Venaille and Sommeria 2007, Du- 
plat and Villermaux 2008), and derived the predictions of these models for the exactly 
unpolluted fraction, P(t), or for the fraction, P(Z c ,t), of the flow with Z below a small 
threshold, Z c . 

To test and constrain the model predictions, we carried out numerical simulations 
evolving decaying scalars in two isothermal turbulent flows with rms Mach numbers of 
0.9 and 6.2. Three passive scalars were included in each flow, and their initial pollutant 
fractions, Pi were set to be 0.5, 0.1 and 0.01, respectively. We found that the mapping 
closure model gives satisfactory predictions for the central part and the high-Z tails of 
the scalar PDF, but underestimates the low-Z tail at large times, especially for scalars 
with small initial pollutant fraction. The left PDF tails become broader with increasing 
flow Mach number, and thus the discrepancy between the mapping closure prediction and 
the simulation results is larger at Mach 6.2. We showed that, in the M = 0.9 flow, the 
scalar PDF is well fit by Gamma distributions at late times, as predicted by Villermaux 
and Duplat (2003). All the closure models adopted in our study were originally developed 
for mixing in incompressible turbulence, and they do not provide successful predictions 
for the scalar PDF in the highly supersonic flow. 

Our simulation results for P(Z c ,t) in the Mach 0.9 flow can be well fit by using 
eqs. (|4.5j) and (|4.9j) from the nonlinear integral model and the continuous convolution 
model of Venaille and Sommeria (2007), respectively. Although these two equations were 
originally derived for the fraction of exactly pollutant-free mass, they provide useful 
fitting functions for P(Z c ,t) with a small finite threshold, Z c . We showed that, for the 
two scalars with Pi 0.1, the evolution of P(Z c ,t) follows the equation P(Z c ,t) — 
P(Z c ,t)\n[P(Z c ,t)]/T con from the continuous convolution model. On the the hand, for 
the scalar with Pi = 0.01, P(Z c ,t) shows different behaviors at early and late times. 
In the early phase, the evolution of P(Z c ,t) is consistent with the equation P(Z c ,t) = 
—P{Z c ,t)[\ — P(Z c ,t)}/Ti n t from the nonlinear integral model, and the later phase is 
well fit by the continuous convolution model. A satisfactory fit to the entire behavior of 
P(Z c ,t) was obtained by connecting the two phases. The continuous convolution model 
starts to apply once the polluted mass fraction is larger than 0.1-0.2. 

When normalized to the flow dynamical time (rd yn ), the decay of P(Z c ,t) is slower 
in the M — 6.2 for two reasons. First, the mixing timescale in units of Td yn , is about 
20% larger than in the M = 0.9 flow. Second, at the same variance, the left tail of the 
scalar PDF is broader at higher Mach numbers. Due to the second effect, the shape of the 
P(Z C , t) curve as a function of time changes as the Mach number increases. We find that 
a generalized version of the self-convolution model ( §3.4 and §4.3; Duplat & Villermaux 
2008) provides a good fitting function, eq. (|4.11j) . to the evolution of P(Z c ,t) in highly 
supersonic turbulence. With n = 4.6, this equation matches our simulation data well for 
the two scalars with Pi 0.1. Like the case of the M = 0.9 flow, we obtained a good fit 
to the simulation result for the scalar with Pi = 0.01 by combining different behaviors 
at early and late times. At early times, we used eq. (|4.5j) from the nonlinear integral 
model, while the later phase was fit by eq. (|4.1ip with n — 4.6. We point out that, 
although it provides a good fitting function to the pristine mass fraction, the generalized 
convolution model does not have a physical connection to how the flow compressibility 
affects turbulent mixing in supersonic flows. Physical PDF closure models are motivated 
to directly incorporate the effects of shocks or the flow compressibility on mixing in 
supersonic turbulence. 
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The fitting functions obtained this study can be used to develop a subgrid model for 
large-scale simulations for mixing of heavy elements in the interstellar media of early 
galaxies. Such a subgrid model would provide an important step toward predicting the 
fraction of pristine gas in the first generation of galaxies. In order to apply our results 
with higher accuracy, we will carry out a systematic numerical study in a future work 
covering a broader range of parameters including varying initial scalar conditions and 
turbulent Mach numbers. 
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Appendix A. The PDF equation 

In this appendix, we derive the equation for the concentration PDF in compressible 
flows using the technique developed by Lundgren (1967). Similar derivations for the 
scalar PDF equation can be found in, e.g., Pope (1976), O'Brien (1980), Dapazo et al. 
(1997) and Pope (2000). Unlike these previous derivations, we adopt a density- weighting 
scheme, which is needed for passive scalar mixing in compressible turbulence. 

Our derivation makes use of a statistical ensemble consisting of many independent 
realizations. We start with the definition of the fine-grained PDF in a single realization. 
Because in a specific realization the concentration field is single- valued at given position 
(x) and time (t), the fine-grained PDF is a delta function 



where Z is the sampling variable. Considering the existence of significant density fluctu- 
ations in supersonic flows, we define a fine-grained PDF with density-weighting 



where the density-weighting factor p is the ratio of the local flow density p(x, t) to the 
average density p. These two fine-grained PDFs are functions of Z, x and t, and their 
dependence on space and time is though C(x, t) and p(x,i). 

We calculate the time-derivatives of q'(Z: x, t) and p'(Z; x, t). Since q'(Z; x, t) depends 
on t only through the quantity Z — C(x, t), we have, 



q'(Z;x,t) = 6[Z-C(x,t)}, 



(Al) 



p'(Z;x,t) = p(x,t)6[Z - C(x,t)}, 



(A 2) 



d g '(Z;x,t) _ dC(x,t) dq'{Z;x,t) 
dt ~~ dt dZ 



(A3) 




(A 4) 



(A 5) 
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We add eqs. (|A 4j) and (|A 5j) . From the continuity equation for p(x, t), the first terms on 
the right hand sides of eqs. JA4| and (|A 5j) add up to zero. Using the advection-diffusion 
equation (|2.1[) of C(x, i) for the sum of the last terms in these two equations, we find, 



dp'(Z;x,t) , d 

- + V-(p(Z;x,t)v) = -^ 



//(Z:x./i ( lv • (p/cVC) + S 



(A 6) 



where we used the fact that, except p'(Z;x, i) or g'(Z;x, i), all the quantities in the 
right-hand-side term are independent of Z. Eq. (|A6|) is essentially a Liouville equation. 
In analogy to the kinetic theory, the concentration here corresponds to the particle mo- 
mentum, and the equation dC jdt = - V- (/O/tVC) + S corresponds to the particle equation 
of motion. 

We now consider the coarse-grained PDF, defined as the ensemble average of the fine- 
grained PDFs over independent realizations. The coarse-grained PDFs with volume- and 
density-weighting are, respectively, g(Z;x,t) — (q' (Z '; x, t)) andp(Z;x,t) = (p'(Z;x,i)), 
where (• ■ •} denotes the ensemble average. The average is over different values of the 
concentration, C(x, t), the flow velocity and density, v(x, t) and p(x.,t), and the scalar 
source, S*(x, i), in different realizations. From eq. (|A 61) . it immediately follows that, 

dp{Z ^ l) + V • <p'(Z; x, t)v) = - A (p' {Z ; x, t) Q V • {pnVC) + sj^ . (A 7) 

The ensemble average terms in eq. (|A 7[) can be expressed in more convenient forms. 
For any stochastic quantity B(x,t), we have the following relation (see, e.g., Pope 2000) 

(q'(Z;x,t)B(x,t)) = q(Z;x,t){B(x,t)\C(x,i) = Z), (A 8) 

where (-B(x, £)|C(x, t) = Z) denotes the ensemble average of B(x,t) conditioned on 
C(x, t) = Z. The conditional mean appears here because the factor q'(Z;x,t) selects 
only the realizations where C(x, t) is equal to Z. Using eq. (IA 8|) . we have 

(p\Z;x,t)B{x,t)) =p(Z;x,t) ^ (x ' f) (A 9) 

where we used p(Z; x, t) = q(Z; x, t)(p\C(x, t) = Z). 

With eqs. (IA 7[> and (I A 9[) . we arrive at the final equation for the coarse-grained PDF 
with density weighting, 

flp(Z;x,t) / (pv\C = Z) \ d f (V-( P kVC)|C = Z) \ 9 / ^gjg = g) 

m + \ p {p\c = z)J dz\ p ( P \c = z) J dz\ p { P \c = z) 

(A 10) 

where we replaced the condition C(x, t) = Z by C — Z for the simplicity of notations. 
Note that the advection term is in a divergence form, which is the motivation for the use 
of a density-weighting scheme in our derivation. A detailed physical discussion of all the 
terms in this equation is given in the text. 
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